Systematic Perturbation Theory for Dynamical Coarse-Graining 
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We demonstrate how the dynamical coarse-graining approach can be systematically extended to 
higher orders in the coupling between system and reservoir. Up to second order in the coupling 
constant we explicitly show that dynamical coarse-graining unconditionally preserves positivity of 
the density matrix - even for bath density matrices that are not in equilibrium and also for time- 
dependent system Hamiltonians. By construction, the approach correctly captures the short-time 
dynamics, i.e., it is suitable to analyze non-Markovian effects. We compare the dynamics with the 
exact solution for highly non-Markovian systems and find a remarkable quality of the coarse-graining 
approach. The extension to higher orders is straightforward but rather tedious. The approach is 
especially useful for bath correlation functions of simple structure and for small system dimensions. 
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I. INTRODUCTION 

The insight that quantum computers may solve cer- 
tain problems such as number factoring [3| and database 
search Q more efficiently than conventional computers 
has given rise to the field of quantum information (for 
an overview see e.g. Q). The conventional paradigm of 
quantum computation relies on unitary operations that 
act on low-dimensional subspaces of the 2"-dimensional 
Hilbertspace of n two-level systems - conventionally 
called qubits. Unfortunately, the dynamics of open quan- 
tum systems is not always unitary such that the im- 
pact of decoherence has to be taken into account. This 
problem also affects alternative schemes such as one-way 
H, 0| , holonomic 0] or adiabatic [1] quantum computa- 
tion. Beyond this, the study of decoherence effects is of 
general interest in the control of quantum systems. 

Often, the dynamics of open quantum systems is an- 
alyzed within the Born-Markov approximation scheme 
[J, 0] . An important criticism raised against this scheme 
is that it does not generallypreserve positivity of the re- 
duced density matrix [l^, HI, EBl , which however is 
necessary for its probability interpretation [iGj . In ad- 
dition, the Born-Markov master equation cannot be ex- 
pected to yield good results for short times, which may 
in the context of quantum computation for example lead 
to false error estimates on required gate operation times 

etc [isini- 

A possible resolution for the latter problem is to study 
non-markovian master equations (that explicitly depend 
on the density matrix at all previous times via a memory 
kernel). However, except for some special cases non- 
markovian master equations are also not guaranteed to 
preserve positivity, and corresponding counterexamples 
can be easily constructed fl^l ■ Technically, master equa- 
tions with memory can for example be solved efficiently 
when the bath correlation functions can be approximated 
by a few decaying exponentials [19]. In the general case, 
they are however difficult if not impossible to solve an- 
alytically. This difficulty transfers to the numerical so- 
lution as well: In order to evolve the density matrix at 



time t, one would generally have to evaluate the solu- 
tion at all previous times t' < t, which corresponds to 
significant computational and storage efforts. 

It is therefore interesting to investigate alternatives 
such as the dynamical coarse-graining (DCG) approach. 
Recently, it been analyzed up to second order in the 
system- reservoir coupling constant (Born approximation) 
j24| . Instead of solving a single quantum master equa- 
tion, the coarse-graining approach defines a continuous 
set of master equations 



(1) 



parametrized by the coarse-graining time r and then in- 
terpolates through the set of solutions at t = r 



Ps{t) = e 



'Pi- 



rn 



Since the Liouville superoperators are of Lindblad 
form [l^l for all t > 0, the second-order dynamical 
coarse-graining approach (DCG2) preserves positivity of 
the density matrix at all times [24]. Note that in the 
general case, the above solution cannot be obtained by 
solving a single Lindblad form master equation merely 
equipped with time-dependent coefficients and should 
therefore be regarded as truly non-Markovian [2l[. The 
conventional Born-Markov-Secular limit is obtained by 
the limit t cx), i.e., ps^^^ = C°°p^^^, whereas in 
the short-time limit, the exact full solution is approxi- 
mated. In addition, it was found for some simple exam- 
ples considered in [2J] that in the weak coupling limit, the 
method approximated the results of the non-Markovian 
master equation for all times remarkably well. 

The purpose of the present paper is two-fold. By in- 
troducing coarse graining in the interaction picture in 
section |TT] we rigorously demonstrate that the method 
will approximate the exact solution for short times by 
construction, i.e., the method is suitable to study non- 
Markovian effects. By including higher orders in the cou- 
pling constant, the agreement between coarse-graining 
and exact solution can be further improved. In addition. 
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we show that up to second order the method uncondi- 
tionaUy preserves positivity of the density matrix, i.e., 
even for bath density matrices that do not commute with 
the bath Hamiltonian or/and for time-dependent system 
Hamiltonians. Wc wiU give several examples for finite- 
size "baths" f subsections IIII Al and IIII B|) . the spin-boson 
model in subsection IIII C[ and we also consider fermionic 
models with transport in subsection IIII Dl 



II. DCG IN THE INTERACTION PICTURE 

A. Preliminaries 

We consider systems where the time-independent 
Hamiltonian can be divided into three parts 



H = Hs + Hi 



SB 



Hp 



(3) 



where denotes the system Hamiltonian, the bath 
(reservoir) Hamiltonian (with [Hs,Hb] = 0), and 



(4) 



couples the two by system (Aa) and bath {Be) operators. 
Note that thereby one has by construction [Aa , B/s] = 
(see section IIIID I for obtaining such a decomposition for 
fermionic systems with transport). 

Note that hermiticity of Hsb = -ffsB imposes some 
constraints on the coupling operators. For example, 
it is always possible to perform a suitable redefini- 
tion of operators by splitting into hermitian and anti- 
hermitian parts (A„ = A^ + A^ and B^ = B^ + B^, 
for system and bath operators, respectively) to obtain 

Hsb = \ {Hsb + ^L) = [A^B^ - ^AtiB^] , such 
that one can always assume hermitian coupling opera- 
tors Aa = Aj, as well as Ba = B^, [4]. For the sake of 
convenience however, we will not assume this form here 
unless stated otherwise. We will use A < 1 as a perturba- 
tion parameter (a-dependent coupling constants can be 
absorbed in the operator definitions). 

In the interaction picture (where we will denote all 
operators by bold symbols) 



+i{Hs + HB)t ^^-f-^^-i{Hs+Hs)t 
Hst 



e • ' "°''p(t)e 
e+'^'^'Aae 



Pit) 

A„(t) 

B„(t) = e+'^«*B„e-*^^* (5) 
the von-Neumann equation reads 

p^^i[HsB{t),pit)] , (6) 
which is formally solved by p(t) = U{t)pQU\t). 

B. Perturbative Expansion 

The time evolution operator in the interaction picture 
is governed hy U — ~iHs-B{t)U {t) , which can be solved 



iteratively. We can define the truncated time evolution 
operator in the interaction picture via 

n t 
fc=0 { 

xSih -h)... Q{tk-i - tk)dh ■■■dtk, (7) 

where the time-ordering is expressed by Heaviside 
step functions. The above operator is unitary up 
to order of A" (assuming that Hsb — ^{X}), i.e., 
Ur^{t)Uj{t) = 1 + e»{A"+i}. Specifically, one has up to 
fourth order 

U^it) = 1 - ^AVl(^) - X'^V2{t) + i\^V3{t) + X'^Viitl^) 
where we can use Eqn. to find for the operators 



Vi(t) 



z 

J2 J dtiA^{h)B^{ti) , 



V2{t) = J2j dtidt2Q{ti-t2)A^{ti)B^{h)x 

a/3 

xAp{t2)Bp{t2), 



I 

V3{t) = Yuj dtidt2dt3Q{ti - t2)Q{t2 



^3) X 



xA^{h)BMAf3{t2)Bf3{t2)A^{t3)B-,{t3), 

t 

Viit) = J2 j ®{ti~t2)Q{t2-H)Q{h-U)x 

xA^{h)Ba,iti)Ap{t2)Bp{t2) X 

xA^{t3)B^{t3)As{U)Bs{U)dtidt2dt3dU ■ (9) 

Using these expressions in the formal solution of the 
density matrix and collecting all terms of the same order 
we obtain 

p{t) = p^-iX[-p^vl{t) + Vi{t)p^ 

+A2 [-poVl(t) + Vi(OPoVl(t) - V2(t)po 
-iX'[p^Vl{t) -V^{t)p,vl{t) 

+ V2{t)p„Vl{t)-V3{t)p, 

+X^ [poVl(t) - Vi(t)poV](i) + V2{t)p,vl{t) 
-V3{t)poVl{t) + V4(t)po] + OiX"} . (10) 

In order to obtain the reduced density matrix, we have to 
perform the trace over the bath degrees of freedom. We 
assume that at io = the density matrix factorizes such 
that we have pg = Tre {po}- Then we can define ps (t) = 
TrB {pit)} and calculate the reduced density matrix at 
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time t 
Ps(t) 



{pgpBvj - Vip°p° + v^pIpIv\ 



+A*TrB 

-V3P°p1?.vJ + V4P°pS^} + 0{A5} 
^ pg + AT,Vg + A2T,V°+A3T3Vg 



(11) 



where we can evaluate the right-hand side by using the 
reservoir correlation functions, see below. 



C. Bath Correlation Functions 

Since we do neither assume a 'priori that the coupling 
operators are hermitian nor that [iJB.Ps] = 0, it IS nec- 
essary to generalize the correlation functions. Denoting 
the index of a hermitian conjugate coupling operator with 
an overbar, we define for the first and second order 

C„(ti) = Tx^{B^{tr)p\) , 
Ca(ti) ^ TTB{Bi{h)pl} , 

a/3(ii,t2) = TrB{S„(ti)S/3(t2K} , 

C^pituh) = TrB{BMBlit2)p°B} , (12) 

and similarly for higher-orders. In terms of these quan- 
tities, the r.h.s. of Eqn. pT|) can be easily evaluated. 



J 



Specifically, we obtain 
t 



^iVg = / dti[-Csiti)p°Al{h) + Co.{h)A^ih)p°] , 



J2 I dhdt2[^C^p{h,t2)eit2-h)p°Al{h)Al{t2) 



+ Gap{tl,t2)Ap{t2)plAl{h) - C„/3(tl,t2)e(ti -i2)A«(il)A^(t2)pg] , 
t 

'^^irf*2rfi3[ + Ca^^(ii,i2,t3)e(t3-i2)e(t2 -ti)pgAl(<i)A^(i2)At(t3) 

a/37 

-Ca^^(tl, t2, t3)e{t2 - ti)A^{t3)p°Al{ti)Al{t2) + Cs,p^{hMM)Q{t2 - t3)Ap{t2)A^{t3)plAl{ti) 

-Ca,f3jih,t2,t3)Qiti - t2)eit2 - t3)A^{ti)Ap{t2)A^{t3)plj , 
t 

J2 J dhdt2dhdU [ + C^-p^-s{h,t2M. U)Q{U - - t2)Q{t2 - ti)plAl{ti)Al{t2)A\{t3)Al{U) 

-C^p^s{tiM. hMWi - t2)e(t2 - h)As{U)plAl{ti)Al{t2)A\{t3) 
+Ca^^5(ii,i2,i3,i4)e(i3 -t4)e(t2-ii)A.^(i3)A5(t4)p°At^(ii)A^(t2) 
-Csf3js{h,t2, h,U)Q{t2 - - U)Ap{t2)A^{t3)As{U)plAl{h) 

+ Ca,(3js{h,t2, h,U)Q{h - t2)e(t2 - i3)e(i3 - U)A^{ti)Ap{t2)A^{t3)As{U)pl 

I 



(13) 



D. Defining the DCG Liouvillian 

It is evident that one can carry on with the expansion 
of the time-evolution operator to arbitrary order in the 
coupling constant A. This will evidently yield good re- 
sults for small A and small times, whereas we would like 
to have a master equation valid for small A and also large 



times. In the original approach [2J| it was shown as a sup- 
portive fact that for f = r the DCG2 solution and the ap- 
proximation pip were equivalent up to OjA^}. Here we 
will demand equivalence between the 7i-th order coarse 
graining solution and Eqn. (jlip at i = r to define our 
perturbation theory. Expanding the Liouvillian superop- 
erator as £^ = A£[ + X^C^ + A^rj + X'^Cl + ©{A^}, 
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we obtain for the solution of Pg(t) = p^{t) at time 
t = T 



Plir) = [l + XrCl 



+ A3 



tCI 



g 1 1 1 



T 

24' 



(14) 



We can clearly match this with equation (fTT|) evaluated 
at i = T order by order to solve for 

1__ . 



Ps 



T 



T 

/"T" 

2 11 



P^ 



g 111 



}pi 



— CTCTCTCT 

2^ 1 1 1 1 



}pg: 



(15) 



where can be extracted from Eqns. 

(fT3)) . Since these equations have to hold for all initial 
conditions Pg we can infer the matrix elements of each 
Liouvillian by comparing coefficients of the matrix ele- 
ments of Pg. 

Equations (jlSp define in combination with (|13p our 
coarse-graining Liouvillian. Evidently, we automatically 
approximate the short-time dynamics of the true solution 
very well by construction with this scheme. 



E. Unconditional Positivity of DCG2 

Here we will show that DCG2 always preserves positiv- 
ity - regardless whether the first order correlation func- 
tions vanish or not. We do not even require [-ffe, Pb] = ^■ 
For simplicity we assume hermitian coupling operators 
Aa = A]^ and = B^,. Then, we obtain from Eqns. 

(Hal 



T 



aP 

-^PsA^{ti)Ap{t2) - ^A^{ti)Ap{t2)ps 

T 

-j^^^y" C'a/3(il,t2)sgn(ti - t2) X 



where we have used Q{x) = i [1 -I- sgn (x)] and 
sgn(— a;) = — sgn(a;) in the last line. From the first of 
the above equations we obtain that the first order Liou- 
villian just generates a unitary evolution 



CIps = -i 



T 

-J2 [ CMAMdh,ps 

^ Q 




.PS 



(17) 



where hermiticity of the Lamb-shift Hamiltonian follows 
directly from hermiticity of the coupling operators (which 
also implies real- valued first order correlation functions). 
In addition, we obtain from consecutive application 



-T-[T-ps] = X J dhdt2Ca^{h)C0{t2) 
X ^Ap{t2)psAaiti) 



^{pS,Ac.{ti)Ai3{t2)} 



This defines the second order Liouvillian as 



(18) 



^2 PS 



■J-^ I Call {ti, t2)sgn {ti -t 



2ri 

xAa{h)A/3{t2)dtldt2, PS 



2 X 



T 

J ^^^1^*2 [Cal3{tl,t2) - C4ti)Cl3{t2) 



af! 



Af3{t2)psAa,{ti) - -{Ao,{ti)Af3(t2),Ps} 



(19) 



The first commutator term induces a unitary evolution 
where hermiticity of the corresponding effective Hamil- 
tonian follows directly from C^^(i2,ii) = Caf3{ti,t2). 
However, in contrast to the standard Born-Markov- 
secular approximation here we have in general 

Hl^,Hs /O. In order to see that the second line 
corresponds to a Lindblad dissipator, we insert identities 
at suitable places 1 = X)a I'') ('^1 obtain 



Q/3 



5 



r;2 
^ah,cd 



ab.cd 



T 

-Y. I ['^o.pihM) - Cc.{h)Cp{t2)] {a\ Ap{t2) \b) (c| \dY dhdh 



(20) 



where we have abbreviated the operators Lab — |a) (^|- The dampening matrix elements can be most conveniently 
evaluated in the energy eigenbasis Hg \a) = Ea\a). 

However, independent from the basis choice it remains to be shown that the dampening matrix is positive semidef- 
inite to get a Lindblad form. In order to see this, we calculate with Eqn. 



E* r,2 
^ablab.cd^cd 



abed 



T 

abed afS q 

T 

^J2Y.<bX^'i J dtidt2Caiti)Ci3it2){a\Ap{t2)\b) (c|A„(ii)M)* 

abed ap q 

r "IF''' 

J2 J B^{ti)x,d (c| \dY dti j Bp{t2)xlk (a| Ap[t2) \b) dt2 

edoc Q abl3 q 

^TrB I j BM^cd (c| A^ih) \d)* dhpl \ttbIyJ Bp{t2)xl^ {a\ Ap{t2) \b) dt2pl 




pI 



i [Tre {K\t)K{t)pI] - Tre {k\t)pI) Tre {K{r)pl]] 
Tre {k\t)K{t)pI] - |TrB {K{t)pI]\ 

I 



(21) 



Whereas the first term in the last line appears positive, 
one might fear that positivity can be spoiled by the exis- 
tence of the additional second term. However, we can 
bound the second term via the Cauchy-Schwarz trace 
inequality [12] |Tr{ylB}|^ < {A^ A] {B^ B] with 

A = K{T)y/^ and B — \ff^ (which exists as is 
positive semidefinite) 



|TrB{if(T)pO}| 



Tre if (r).MJpO 



< Trp 



IpIk\t)K{t). 



'pUpb 



xTre 

= TrB{KHT)Kir)p'^} 



(22) 



Remembering that Tre {K^t)K{t)p'^} > for any op- 
erator K{t) we therefore obtain for r > 



Not only is the case Ca{ti) ^ considered but in ad- 
dition, we do not constrain ourselves to bath density 
matrices in thermal equilibrium, i.e., one also has pos- 
itivity for [pq , Hb] 7^ . It is an interesting avenue of 
further research to compare DCG with other methods 
within the context of nonequilibrium environments ^2^. 
Beyond that, all of the above arguments go through if 
the system Hamiltonian is time-dependent. In this case, 
the coupling operators in the interaction picture have to 
obey Aa — +i [Hs{t), Aa{t)], such that the challenge is 
then to calculate the matrix elements (a| Aot{t) \b). 

Under the assumptions Ca(t) — (no first order 
correlation functions), C'ap(ti,t2) — Capiti — t2) = 
TiB {Ba{ti — t2)BppQ^ (reservoir in thermal equilib- 
rium) we can insert the Fourier transforms of Cap{ti —t2) 
and Cajsiti — t2)sgn {ti — 12). If in addition the system 
Hamiltonian is time- independent, we may calculate the 
time integrals analytically. Then, we can make use of the 
identity for discrete a, b (see e.g. appendix F of ref. [24|) 



Y <blab,cd^cd > , 



(23) 



abed 



lim rsinc 

r — ►oo 



{iO + a)T 



{ijj + b)T 



2TrSabSiuj + c(j2A) 



i.e., we have generated a Lindblad form master equa- 
tion. This result goes beyond ref [24^] in several aspects: 



to calculate the large time limit of the DCG2 approach. 
In complete analogy to ref [IJ] we obtain the Born- 
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Markov secular approximation [4| in this limit. 

Unfortunately, the unconditional preservation of posi- 
tivity is not preserved by higher orders within DCG (al- 
though of course, in the weak coupling limit the nice 
properties of DCG2 will dominate). 



III. EXAMPLES 

In the following, we will test the DCG approach with 
simple examples for which at least in special cases an an- 
alytical solution exists. For finite-size reservoirs the cor- 
relation functions are non-decaying and these systems are 
inherently non-Markovian (exhibiting for example recur- 
rences) , cf. the examples in subsections IIII Al and IIIIBI 
For quasi-continuous reservoirs we will compare the per- 
formance of the DCG approach with the Born-Markov 
approximation, see subsections IIII CI and IIII Dl 



A. DCG2 for two spins 

We consider a highly non-markovian system (S) by us- 
ing a very small reservoir (B), namely just a single further 
spin 



B J 



Hs = 

HsB = A(Ts • CTB = A [cTg cr| -I- cr^ ® cr]^ -I- (Tg 



(25) 



i.e., the index of the coupling operators runs from one 
to three. Note that all coupling operators are hermi- 
tian, such that we may omit overbars and daggers in 
Eqn. ([131). We assume that the initial bath density 
matrix is diagonal in order to simplify all expressions 



- 1 1 - pj^O 



The exact solution can be ob- 



tained by exponentiating the Hamiltonian and tracing 
out the bath spin (not shown). As in subsection (jll E[) we 
decompose the Liouville operator into unitary and non- 
unitary contributions, where we have first and second 
order contributions in the unitary action of decoherence 
iJgjj — + and second order contributions for 

the dissipative action Yab,cd = lab,cd- 

Transforming the coupling operators into the interac- 
tion picture we obtain Bi{t) = cos(2f2t)crg — sin(2f7t)(Tg, 
Ax{t) = cos(2wt)crg - sm{2Ljt)al, B-2{t) = cos(2m)(T^ -|- 
sin(2m)a|, A2(t) = cos(2wt)cr^ sin(2wt)CTg, B^it) = 
CTg, and As{t) ~ o-g. From this, we obtain the time- 
independent first order correlation functions 

Ci{t) = , C2(t) = , C^{t) = 2pS^° - 1 , (26) 

which yield for the first order Lamb shift Hamiltonian 
from Eqn. ([T71) 



H. 



-r,l 



eff 



A(2p°" - l)<j, 



s • 



(27) 



The non-vanishing second order correlation functions 
equate to 



Cn 

Cl2 
C33 



cos[2(ii - t2)^] - i(l - 2^5^") sin[2(ti - 12)^] , 

- 2^^") cos[2(ti - ^2)^^] - sin[2(ti - ^2)^^] , 
i{l - 2p°^) cos[2(ii - t2)n] -f sin[2(ii - 12)^] , 
cos[2(ti - t2)^] - i(l - 2^5^*^) sin[2(ti - 12)^] , 
1, (28) 



where we have omitted the time-dependencies for brevity. 
This can be inserted in the expression for the second 
order Lamb-shift Hamiltonian in Eqn. to yield 



H 



eff 



2A2 



{1-sinc [2T(f7-cj)]} X 



(29) 



which commutes with the system Hamiltonian. 

The second order dissipative terms must be 
calculated from the dissipative parts of Eqn. 
(|19p . where we obtain for the non- vanishing 
matrix elements of the dampening matrix 



T,2 



loom - 4AV(1 - pM°, loo.n = -4AV(1 - p^")pO°, 
7itoo = -4AV(1 - pO")pO", 7itii = 4AV(1 - p^'. 
7^|2^^ = 4A2rp00sinc2 [r(f} - lu)] , and 
7i6^io ~ 4A^r (1 — pg*) sinc^ [T(il — w)], which shows 
(e.g., by the Gershgorin circle theorem [^) that ^ab,cd is 
positive semidefinite. The solution of the coarse-graining 
master equation Pg(i) = p^{t) can be conveniently 
obtained by exploiting that diagonal and off-diagonal 
matrix elements decouple. From the diagonal equations 



7oi,oiPs (^) - 7i6,ioPs W 



-r,2 
7oi,01 



T,2 . T,2 
7oi,01 + 7lO,10 



4AVpg'sinc'' [T{n - Lu)] 
-4AVsinc2 [T{n-uj)] p™(t) 



(30) 



we obtain the solution at r = i 



IpIWIo) 



exp < —4 



A2 



■sin^ m^u)]\p's"iO) 



4A2 sin^ [t{n - Lo)] 



1 — exp < — 



PB : 



which does admit for complete recurrences of the popu- 
lations, see figure [T] (a). 
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For the ofF-diagonal equation 



{ 



T,2 _ 1 / r,2 , r,2 . t,2 , t,2 \ 

7oo,ii 2 l^'^oo.oo + 7oi,oi + 7ioao + 7ii,ii j 

(ii h:^' ii) + {i\h:^^ ii) - (oi h:^' io) 

-(0|fZ-:^=^|0>)]p°Ht) 

8AVpS^"(1 - p°^) - 2AVsinc2 [t(17 - u)] 
,/2A2(l-sinc [2T(f7-tj)]) 



+ 



Q — uj 
2A(l-2p^°))}p°H0 



(31) 



we obtain the solution 



Pj("-)l) +2A(l-2p°°)j 

e L J pui 



where we will generally observe a decay whenever 
pOO) ^ 0, see figure [11(b). 



B. DCG4 for two spins 



In order to keep the calculations for DCG4 very simple, 
we consider 



Hs - ujal , Hb = nal, , Hsb = Aag , (33) 

where also here the coupling operators are hermitian. In 
this case, the bath correlation functions are all time- 
independent, which enables a convenient calculation of 
the Liouvillian matrix elements. The example is of course 
a bit trivial, since the exact solution for the reduced den- 
sity matrix does not depend on $7. Note however, that 
unlike the pure dephasing limit considered in [26| this 
case still holds some time-dependence that can be found 
in the system operator in the interaction picture. 

The exact solution can be calculated by exponentiat- 
ing the complete Hamiltonian and tracing out the second 
spin in the solution for the density matrix as in sub- 
section nil Al and in a similar manner we determine the 
DCGl, DCG2, DCG3, and DCG4 solutions by directly 
determining the 4x4 Liouvillian matrix as described in 
subsection III Dl (not shown). The solution is then ob- 
tained by exponentiating the Liouvillian. The resulting 
solution for the diagonals is displayed in figure O (a) and 
for the off-diagonals in figure [5] (b). 




time t [a.u.j 

FIG. 1: [Color Online] Comparison of exact (solid black) and 
DCG2 (dashed red) solutions for the diagonal (a) and off- 
diagonal (b) matrix elements of the density matrix. Imagi- 
nary parts are displayed with thin lines. Naturally, the exact 
solution displays complete recurrences. For the diagonal ma- 
trix elements, this feature is well reproduced by the DCG2 
approach, whereas for the off-diagonals only the short-time 
dynamics is well approximated. The other parameters have 
been chosen as follows A — 0.25, ui = 1.0, Q = 2.0, and 
= 0.5. 



C. Spin-Boson model 

We consider a single system spin coupled to a bath of 
bosonic modes {ujk > 0) 



SB 



AA. 



hkbk + hlb\ 



(34) 



where for simplicity we have restricted ourselves to the 
case of single-operator coupling, and the operator A = A'^ 
will be specified later-on. 

We will consider a thermalized initial bath density ma- 
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FIG. 2: [Color Online] Comparison of exact (solid black) and 
DCGl (dotted orange), DCG2 (dashed red), DCG3 (long- 
dashed green), and DCG4 solutions (dot-dashed lines) for the 
diagonal (a) and off-diagonal (b) matrix elements of the den- 
sity matrix. For small times, DCG4 is superior to the coarse- 
graining methods of smaller accuracy. By construction, the 
exact solution shows complete recurrences. For larger times, 
all coarse graining methods also display recurrences but all of 
them miss the exact solution (not shown). Parameters have 
been chosen as A = 0.5, a; = 1, and pg" = 1. 



trix 



TrB{e 



-/3i?E 



} 



(35) 



where (3 = (fceT) ^ denotes the inverse reservoir temper- 
ature. 



must be balanced for all modes to obtain a nonvanish- 
ing result. Therefore, we conclude (since only one op- 
erator is involved, we may omit the indices) C{t\) = 
= C{t\,t2,t-i). In the interaction picture, the anni- 
hilation and creation operators transform according to 
hk(t) = Q+'i.Hst^^^-iHBt ^ e-'^^'bfc and the hermitian 
conjugate, respectively. 

The second order correlation function then evaluates 

to 



^+iLo{ti—t2) 



— j dLoG[ijj)^n[Lo)e 



+ [l + n(u;)]e-*"(*i-*^'} 

(36) 



4-00 

l_ f G{\lo\) 
2tt 



where the bosonic occupation number is given 
by n{u!) = ■ In the above equation, we 

have assumed a quasi-continuous spectral density 
G{uj) = 2TrJ2k\^k\^^{^ — ^k) to convert the sum into 
an integral. When we parametrize the spectral density 
as 



(37) 



where lUc denotes a cutoff frequency and the parameter 
S governs the slope at a; = 0, we can obtain an analytic 
solution for the correlation function [23, [28l 



C(tl,t2) 



Gor(l + 5) 



+C (1 + ^,1 



C 1 + "5^:75—+* ^ 



1 _ . ftl -<2) 



(38) 



in terms of generalized Riemann Zeta functions C,{x,y). 

The next non-vanishing correlation function is fourth 
order, where we obtain 

C{h,t2,t3,U) = C{t2,t3)G{ti,U) + C{h,t3)C{t2,U) 

+G{h,t2)G{t3,U), (39) 

which can for example be obtained usin g W icks theorem 
(for a special case see also Eqn. (61) in jll|). 



1. Bath Correlation Functions 

We evaluate the traces in the correlation functions in 
Eqn. (fT2)) in Fock-space, where the bath density ma- 
trix ([55]) is diagonal. By doing so it becomes obvi- 
ous that the number of creation and annihilation opera- 
tors in each term of the bath-bath correlation functions 



2. Pure Dephasing 

The case of pure dephasing A = cr^ is exactly solv- 
able [1^ [13] and it is known that DCG2 already yields 
the exact result |23|. The exact solution predicts time- 
independent diagonal matrix elements and a decay of the 
off-diagonal matrix element according to (in the inter- 
action picture, cf. Eqn. (82) in 30] in the limit of a 
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continuous bath spectrum) 



we obtain for the fourth order contribution 



Poi{t) 

m 



e-"(*Voi(0), 



8A^ 
2tt 



^, ,sin2(wt/2) , 
G(w) \ coth 



duj ,(40) 



where the additional factor of ^ in comparison to [l^l 
results from a different definition of the spectral density 
G(aj). By taking the time derivative of the exact solution 
density matrix we obtain a closed master equation that is 
not of Lindblad form (not even with time-dependent coef- 
ficients) but nevertheless must - as it is exact - preserve 
positivity. Hence, one would regard this case as truly 
non-Markovian [21|. The corresponding steady state is 
also derived by the equation of motion method in ap- 
pendix [X] For pure dephasing we have A(t) = A — a^. 
We observe a decoupled evolution of diagonal and off- 
diagonal matrix elements of the density matrix. 

For the second order contribution we have (using 
Oih - t2) + Q{t2 - h) = 1) from Eqn. ^ 



T 

^2^8 = / C{hM) yplo' - pI] dhdt2 , (41) 



from which we obtain 



(0|T-pg|0) 
(0|T/p°|l) 



(l|T2-pg|l)=0, 

r 

-2 j C{tiM)dtidt2pf 



oo 

8 /■ . . sin^(W2) 



coth 



which leads to the same exponential decay as with the 
exact solution ([3D]), i.e., as noted earlier [2^, DCG2 yields 
the exact solution in this case. Note that this example 
demonstrates explicitly that the DCG2 solution cannot 
be generally written as the solution of a single Lindblad- 
type master equation with time-dependent coefficients, 
such that it should be classified as non-Markovian [2l| . 

Using the relation 



(0|T/p°|0) = (l|T/p°|l)=0, 

r 

(0|T/pg|l) = 2pl^ j dtidt2dhdUC{ti,t2,h,U) X 



Q{h~t2)e{t2-ti) 



= 2 / dtidt2dt:idUC{tiM)C{HM)PW^) 



where we have exploited the symmetries of the fourth 
order correlation functions under exchange of the argu- 
ments and the relation 



2 = 



+ Q[t2 




-U) 


+ e{t3 


-t2)Q{t2 


-h) 


fe(ii - 


t2Mt2 ~ 


u) + 


Q(t2- 


ti)eih - 




f e(i3 - 


t2Mt2 - 


u) + 


e{t2- 


i3)e(t3 - 


h) 


f e(t4 - 


hMh - 


t2) + 


eih- 


u)e{u - 




f e(i3 - 




t2) + 


e{u- 


^3)6(^3 - 


h) 


fe(ti - 


U)Q{U - 


t2) + 




h)eih - 


h) 



The non-vanishing off-diagonal contribution has to be 
compared with the counter-term arising from the second 
order 



i(0|T,nT,-pg]|l) = i 



-2 j C{ti,t2)dtidt2 



C{tiM)dtidt2 



Ps 



(7(^3, t^jdt^dt^ 



pgi (45) 



Since the diagonal elements of the density matrix are 
neither affected by nor by we conclude that we 
have for pure dephasing 



-'4 ^ 2 2 -'2 ' 



(46) 



such that DCG4 will yield the same result as DCG2. 
Since DCG2 is already exact for this case, this cancella- 
tion is a strong indicator for the correctness of our fourth 
order correlation function ([5^ . 



= +eih ~ h)e{t3 ~ t2)e{t2 ~ h) 
+e{t3 ~ h)e{t2 ~ h) 
+e{h-t2)e{t2-t3)e{t3-ti) 

-e{t2 - hMh - U) - Q{t3 - t2)Q{t2 ~ h) (42) 



3. Dissipation 

We are now in a position to apply DCG to more inter- 
esting coupling operators (picking up a time-dependence 



10 



in the interaction picture) that also affect the evolution A{t) = cos{edt)a^ + sin(edi)o'^. Inserting this result we 
of the diagonal elements of the density matrix. Trans- find for the second order term in Eqn. (jl3p 
forming A ~ into the interaction picture we obtain 



I^2"[PS] |0) = / dhdt2C{ti,t2 



g-ied(ti-t2)p00 _|_ g+ied(ti-t2)pll 



r 



and similarly for the other terms, such that by arranging the matrix elements of the 2x2 density matrix in a 4- 
dimensional vector as (Ps"' Ps^' Ps"' Ps^) '^^'^ matrix elements of the corresponding 4x4 superoperator to 
be 



I _g-»ed(ti-t2) Q Q _j_g+jed(ii-t2) \ 

Q _g-ied|tl— 12| _j_g — i£d(tl+t2) Q 

Q ^g+i£d(*l+t2) _g+ied|tl-*2| Q 

y _|_g-ied(ti-t2) Q Q _g+ied(ti-t2) 



dt\dt2 



(48) 



r 



such that we observe a decoupled evolution of diagonal which corresponds to the thermalized system density ma- 



and off-diagonal matrix elements. Defining 



mii(T) 



mi4(r) 



77141 (t) 
77144 (t) 



' G(|u;|) 



r 
2^ 



~oo 



2 /14 
2 

T 
'2^ 



|g/3^ - 1| 

> 

|g/3a. _ 1| 



smc 



(w - Ed): 



(w + £d)i 



du} . 



dw . 



(^2")41 = -"^llW: 
(^2^)44 = -"^14(t") 



(49) 



we obtain the second order solution for the diagonals (us- 
ing trace conservation) 



POO(^) = Poo CXP [A^ (mil(T) - 77il4(T))] 

1 - cxp [A^ (mii(T) - mi4(r))] 



Y _ ™il(r) 
mi4(r) 



(50) 



trix that consistent with our expectations, compare also 
appendix [XI The same stationary state can also be ob- 
tained by using the method of equation of motion and 
truncating correlations at second order between system 
and reservoir. 

The Markovian limit is obtained by using 



m^^{t) = t hm -7nn(T) = -tG(|ed|)|?^(+ed)| , 

r— >-oo T 

m^f (i) = t hm -mi4(r) = +iG(|£dl)k(-£d)| , (52) 



where we have again used identity ([M)) in Eqn. (|49p . 
Evidently, the above equation leads to the same steady 
state as the Born-Markov approximation ([5T|) . 



By virtue of Eqn. ([1^ and using that 

[C(ti, i2, ^3, ^4) = ^(^4, is, i2, ii)]* we obtain for the 
diagonal part 



In Eqn. (gH) it is already obvious that for finite times (Oj T^'^Pg |0) = 2 / (c(ti, ^2, ^3, i4)e~*^'*^*^~*^+*^^*'*^ | 
T all frequencies will contribute to the matrix elements ^ 
77iij(r) in contrast to the Markov approximation, where 
only G'(£d) is relevant. With using that the bandfilter 
sine functions transform into Dirac-Delta functions in 
Eqn. we can perform the limit r — » 00 to obtain 

the steady state 



6(^2 - ^3)6(^3 - t4)dtidt2dt3dt4p°° 

T 

-2 / sft|c(ii,i2,i3,i4)e+*"<^(*^~*'+*-'"*'') 



1 



Poo 



1 -I- e-^=d 



(51) 



6(^2 - ^3)6(^3 - t4)dtidt2dt3dt4pl^ 
Pn(r)p°°+Pi4(T)pii. 
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This result has to be combined with the counterterm aris- 
ing from the squared second order contribution. Defining 

A** 

Pii(t) = A^mii(T) - Y + ™i4(t)to4i(t)] 

+AVi(t), 

A"* 

Puir) = A^mi4(T) - Y [™ii('^)™14(t) + toi4(t)to44(t)] 

+AV4(t), (54) 

we therefore obtain the fourth order solution 

PooM = Pooexp[pii(T) -pi4(t)] 

1 -exp[pii(T) -Pi4(t)] 



1 



£ll(l) 
Pl4(-r) 



(55) 



A general exact solution is unfortunately not available 
for this case. However, it is interesting to note that when 
considering the Markov limit /3 = 0, Wc — > oo, and 5 = 
(where the Markov approximation becomes exact) the 
correlation function (|36p becomes a (5-function and we see 
that in this limit, the fourth order term is cancelled by the 
squared second order counter term. For comparison, we 
plot Born-Markov solution, DCG2, and DCG4 solutions 
in figure [3] 

For the dissipative spin-boson model and Ohmic dissi- 
pation (S — I) one obtains an exponential decay of the 
expectation value (a^) (t) in the long time limit [3l| (note 
the rotation and ct^ "O"^)- This corresponds 

to a decay of the off-diagonal matrix elements and is of 
course also reproduced by the DCG approach. 



D. Fano- Anderson Model 

We consider the Fano- Anderson model [H, [s^: two 
leads that are connected by a single quantum dot, 
through which electrons may tunnel from one lead to the 
other. The Hamiltonian is given by 

H = Hs + Hb + HsB , 



HsB 



ka 



(56) 



with ferniionic operators creating an alectron with mo- 
mentum k in the left or right lead a G {L,R} for c^^ or 
in the quantum dot d^. Due to the fermionic anticom- 
mutation relations, the model (|56|) can be solved exactly 
for example with Greens functions [34, 3§|. We provide 
a simplified derivation based on the equation of motion 
method in appendix [Bl 

In contrast to our assumptions in section [TTl the op- 
erators d and Cka do not act on separate Hilbert spaces, 
which is evident from their anticommutation relations 



DCG2 
BMS 



/ y/ 



1 10 

time t [a.u.] 



8 o-' 

Q. 



DCG2 
DCG4 
BMS 



\ . 



w. 



time [a.u.] 

FIG. 3: [Color Online] Comparison of DCG2 (dashed red), 
DCG4 (dot-dashed blue), and BMS (dotted green) approxi- 
mations to the dissipative spin-boson model. In figure (a) we 
consider different temperatures /3 = 0.2ed (bold), /3 = LOsd 
(medium), and /3 = S.Osd (thin) in the weak coupling limit 
=0.1 and we see that the steady states always coincide 
whereas for small times DCG and BMS solutions differ. In 
figure (b) we only show the short time-dynamics for A^ =0.1 
and /3 = l.Osd- The other parameters have been chosen as 
£d = 1, a;c = 1, Go = 1, 5 = 1. 



|d, cj[,jj| = 0. This becomes even more explicit via the 
decomposition 



d = |0)(1|®1, 

Cka = [|0)(0|-|l)(l|]®Cfca, 



(57) 



where |0) and |1) denote the empty and filled dot 
states, respectively, and the fermionic operators Cka act 
only on the (distinct) Fock space of the remaining sites 

in the leads with |cfco,4'a'} = 4fe"^aa'- Naturally, 
the above decomposition obeys the original anticom- 
mutation relations such as \d^Cka\ = and we con- 
clude for the operators in the interaction Hamiltonian 

C^cL = - |0>(1| ® 4a' and dtcfea = + |1) (0| ® Cka, SUCh 

that now the new operators commute by construction. 
Similiar decompositions are also possible systems con- 
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taining more than one site. We identify 

A, = -|0)(1| , A2 = -|l)(0| , 



hermitian in this case. The correlation functions relevant 
for the evolution of the diagonal matrix elements can be 
calculated by introducing continuum tunneling rates via 



Bl = Y^^kacla^ B2=E^feaCfea. (58) H = 2^ ^fc l^^a | - ^.a ) 



ka ki 

Wc assume that there are no correlations between left 



+ 00 



vve assume Liiai mere are no correlations ueiweeii leiL p 

and right leads. Here we will just consider the infinite Ch{ti — ^2) = Ci2{ti,t2) ~ 2^ rL(w)e^*'^''*^ ^"^^du). 
bias limit (although this is not crucial, it enables for an J^^ 



analytic calculation of all integrals) : Taking the chemical 
potentials to plus or minus infinity for the left and right 
leads, respectively (/j,l +00 and /iR —00), we obtain 
for the fermionic occupation number 



1 



/ g/3(lJfc-/JL) ^ 1 



1 



1. 



0. 



(59) 



We observe that the coupling operators are non 



+00 

C^{t2~h) = C2i{ti,t2) = ^ J rR(a;)e-*'^(*i-*^)dc^, 

—00 

(^1212(^1,^2,^3,^4) = Ch{ti — t2)CL{t3 — t4) 

+Ci^{h-u)CK{h~h), 

^2121(^1,^2,^3,^4) = C'r(^2 ^ il)C'R(^4 ^ ^3) 

+CK{u-ti)c^{t2-him) 



J 



From Eqn. p3p . we can derive the second order approximation 



^2"PS = j C?ilrf^2[-C2l(tl,t2)e(t2-il)pgAl(ti)A|(t2)-Ci2(tl,t2)e(i2-il)pgA^(ii)A{(t2) 


+C2i{h,t2)A^{t2)plA\{h) + Ci2{tut2)A-z{t2)plAl{ti) 

~Ci2{h,t2)Q{h - t2)Ai{ti)A2it2)pl ~ C2l{tl,t2)Q{tl - ^2 ) ^2 (tl ) Ai (i2)p° 

+ 00 r 

= / ^rR(c.) j diidt2e-'(-^^»*-*^) [-pI |1) (1| 6(^2 - h) + |0) (1| pI |1) (0| - |1) (1| plQ{t, - t2)] 

-00 

+OC r 

+ / ^rL(c^) J diidi2e+^(--^^)(*^-*^) [-p° |o) (o| e(t2 - h) + |i) (o| p° |o) (i| - |o) (o| p°e{h - 12)^61) 

-00 

from which we can infer the matrix elements of the second-order Liouvillian viaC^pl = ^T^p%- The Born-Markov- 
secular approximation is obtained by the Liouvillian £^ with the help of Eqn . 

I 

When we parametrize the tunneling rates by ments of that govern the evolution of the diagonals 
Lorentzians 13611 



^0 4t/j- \ /it/ 



, rL(^) = 



, (62) 



/mn(r) m^jr) \ ^ f dt,dt2e-^^^^'^~'^^ x 
y TO4i(t) m44(T) J J 



we obtain an analytic expression for the bath correlation 
functions in terms of a single decaying exponential 



'Ci,{h-t2) +Cn{ti-t2) 

+ CL{tl~t2) -C^{tl-t2) 



Then, we obtain from (Ojl^'^pg |0) for the evolution of 
the diagonal matrix element 



CR(t) = I^e-I*l^''+'^^*, CL(t) = ^^e-l*l^^+^^^*, (63) Poo - 



™ii(t) / ,n , "^i4(t) . ■ 



(65) 



where we can exploit trace conservation 
such that we can analytically calculate the matrix ele- p\^ (<) = 1 — Pgg [t] (this feature is trivially fulfilled 
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by the coarse-graining approach). Afterwards, the above 
equation can expUcitly be solved for 

Poo(*) = Poo(0)cxp [mii(T) - toi4(t)] ^ 
1 - exp {A^ [mii(T) - mi4(T)] ^} 



2 _ ™ii(t) 
m-nij) 



(66) 



For the fourth-order contribution we obtain (extensively using (^1212(^1,^2,^3,^4) = C'2i2i(ii, ^2, is, ^4), 
Ciii2{ti,t2,h,ti) = etc.) 

r 

^/PS = j dildi2dt3rfi4[+C212l(il,<2,t3,<4)e(t4-i3)e(i3-t2)e(t2-il)pgAl(ti)A|(t2)Al(t3)A|(t4) 


+Ci2i2(ii, i2, h,U)Q{U - t3)e(t3 - i2)e(i2 - h)plAl{ti)A{{t2)Al{h)Al{U) 
-C2i2i{ti,t2,h,U)Q{h - t2)Q{t2 - h)A^{U)plA{{h)Al{t2)A\{t^) 
-Ci2i2{ti,t2,h,U)Q{h - t2)e(t2 - ti)A^{U)plAl{h)A{{t2)Al{t^) 
+C2ii2{ti,t2,h,U)Q{h - U)Q{t2 - h)A^{h)A^{U)plA{{h)Al(t2) 
+C2i2i{ti,t2,h,U)Q{h ^ U)Q{t2 - ti)A^{h)A^{U)plA{{h)Al{t2) 
+Ci2i2{ti,t2,h,U)Q{h ^ U)Q{t2 - h)A^{h)A^{U)plAl{h)A\(t2) 
+Ci22i{ti,t2,h,U)Q{h - U)Q{t2 - h)A^{h)A^{U)plAl{h)A\(t2) 
-C2i2i{ti,t2,h,U)Q{t2 - - U)A^{t2)A^{h)A^{tA)plA{(ti) 

-Ci2i2{ti,t2,h,U)Q{t2 ~ - U)A^{t2)A^{h)A^{U)plAl(ti) 

+Ci212{h,t2,t3,U)eih-t2)eit2~t3)eit3~U)Aiiti)A2it2)Ai{t3)A2iti)p° 

+ C212l{h,t2,t3,h)e{ti - t2)eit2 - t3)e{t3 - U)A2{h)Ai{t2)A2it3)Ai{U)p°\ . (67) 

The relevant part in above equation for the evolution of the diagonals evaluates by using Eqn. (j42|) to 
(^4"Ps)ii = PiiMpSoW+Pi4(t)p[i(<), 

r 

pii(r) = + / dtidt2dt3dUe~'^''''+*^-*^+*^-*''^ x 







[CLih - <2)Cl(<3 - h) + Ci^ih - h)Cn{t3 - 12)] [e(t3 - t2)e(i2 - h) + e(i2 - t3)e(t3 - ^4)] 

r 

Pi4(t) = - J dtidt2dt3dUe+''^^''-''+*'-'^'> X 



[Cn{t2 - h)Cn{t4 - ts) + Cr{U - ti)Ci^lt2 - is)] [e(i3 - ^2)6(^2 - h) + Q{t2 - t3)Q{t3 - U)] , (68) 
where again all integrals can be solved analytically yielding even lengthier expressions than before. 

Together with we obtain from 



\2 \4 / 2 



(0| T^T^pI |0) 



[rnii(r)rnii(T) 

[TOll(r)TOl4(T) 



mi4(T)m4i(r)] 

TOl4(r)TO44(T)] 



1 



= - A^T^- + A^T/ - -T/T^- 69) 
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the following differential equation for the diagonals 

A* 



P6o 



1 r A r 

A^toii(t) - — TOii(r)TOii(r) 

T Z I 

+mu{T)mii{T) + A'*Pii(t)|pJo(/:) 



+ - A^mi4(T) 



mii(T)mi4(T) 



+mi4(r)m44(T)] + X^puir)} pl^{t) . (70) 

The above equation can be explicitly solved for pJq(t) 
by imposing trace conservation as in subsection IIII CI 
It can be shown analytically that in the flatband limit 
(5r — > oo, (5l — > oo (where for infinite bias the Markovian 
approximation becomes exact), the fourth order correc- 
tion in the above equation is cancelled by the counterterm 
from the second order for all graining times r. Moreover, 
one can also show analytically that the apparent diver- 
gence for large graining times pn oc and pi4 cx 
is precisely cancelled by the fourth order counter terms 
arising from the second order. 

In contrast, one obtains under the Born-Markov (the 
secular approximation has no effect for this particular 
example) approximation the solution 



Poo(t) = 



rR(ed) 



1 



,-A2[rL(ed)-i-rR(£d)]-> 



rL(ed)+rR(ed) 
+Poo(0)e-^'F-(^-)+r-(--)l^ 



(71) 



which has been derived using Eqn. (|52p and identity 
p4)) . From the Lorentzian tunneling rates ((62)) we see 
that the Markovian solution is completely independent 
on the width of the tunneling rates (5r , 5l, see also fig- 
ure m Especially for small widths 5l the correla- 
tion functions (|63p decay very slowly and we also ob- 
serve a large difference between Born-Markov and exact 
solution, whereas the DCG solutions perform comparably 
well. With the exact solution from appendix [B] wc plot 
the Born-Markov solution, and DCG2 as well as DCG4 
solutions in figurelHlfor varying coupling strengths as well 
as for different model symmetries. 



IV. SUMMARY 

The dynamical coarse-graining approach has been ex- 
tended to higher orders. By performing the derivation 
in the interaction picture, we have directly demonstrated 
that the DCG solution must approximate the exact so- 
lution by construction for small times. This has been 
confirmed by several examples. Interestingly, the DCG 
method even reproduced the complete recurrences of the 
diagonal density matrix elements in case of small reser- 
voirs. For continuous reservoirs, the short time dynamics 
of DCG4 has always been superior to the short-time per- 
formance of DCG2. This however need not be the case for 
the large time limit. However, the performance of DCG2 
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■ — ■ DCG4 

■ ■ ■ ■ Born-Markov 
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FIG. 4: [Color Online] Comparison of exact (solid black), 
DCG2 (dashed red), DCG4 (dot-dashed blue), and Born- 
Markov (dotted green, same for all parameters) solutions for 
the Fano- Anderson model. The bold lines (Jl = 5r = 0.1) 
show the highly non-Markovian regime whereas the medium 
thickness lines ((5l = 5r = 1.0) denote a regime where the 
Markovian approximation performs comparably well. The 
thin lines (Jr = 25l — 0.1) demonstrate the failure of the 
DCG solutions in the large-time limit. The other param- 
eters have been chosen as = 0.1, F^ = Fl = 1, and 

Sr ~ £d ~ £h ~ 1. 



is always better (short time) or equal (large time) to the 
performance of the Born-Markov-secular approximation. 

We have shown that DCG2 unconditionally preserves 
positivity. Going beyond [24] , this also includes reservoirs 
that are not in equlibrium. In addition, the presentation 
in the interaction picture leads to a much simpler form 
of DCG2, such that now it appears at least as simple (if 
not simpler) than the conventional Born-Markov theory. 
Positivity is not unconditionally preserved for higher or- 
ders of DCG. 

Unfortunately, the DCG method is computationally 
quite demanding in the interesting case of large (con- 
tinuous) reservoirs, as it requires the evaluation of high- 
dimensional integrals. As the dimension of some integrals 
can be reduced by analytical integration, the efficiency 
of DCG may therefore strongly depend on the structure 
of the bath correlation functions. For systems that are 
larger than the single spins considered here, it will also 
prove difficult to calculate the exponential of C^t for all 
t and T with moderate computational effort. If however 
the result is only of interest at a specific time (as is the 
case for example in adiabatic quantum computation) one 
only has to exponentiate a single matrix or - even simpler 
- evolve the density matrix according to a single Liou- 
villian. The fact that DCG2 unconditionally preserves 
positivity also for time-dependent system Hamiltonians 
renders the method a good candidate for analyzing the 
corrections of decoherence to adiabatic quantum compu- 
tation ^ without the necessity of reverting to conven- 
tional Born-Markov-secular theory [sj ■ 
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time t [a.u.] 



FIG. 5: [Color Online] Comparison of exact (solid black), 
DCG2 (dashed red), DCC4 (dot-dashed blue), and Born- 
Markov (dotted green) solutions for the Fano- Anderson 
model. In figure (a), we have considered symmetric maxi- 
mum tunneling rates Fl ~ Fp = 1 for the weak coupling limit 
= 0.1 (bold lines) and the strong coupling limit = 1.0 
(thin lines). It is visible that the steady state of BMS and 
DCG2 solutions does not depend on the coupling strength 
and in the strong coupling limit, the steady state of DCC4 
might actually be worse than that of DCG2. In figure (b), we 
have used A^ = 0.1, F^ — 1.0 and varied the left maximum 
tunneling rate as Fl ~ 2.0 (bold), Fl = 5.0 (medium), and 
Fl = 10.0 (thin). For small times, DCG4 always yields a bet- 
ter result than DCG2 and both DCG solutions are better than 
the BMS approximation. The latter performs particularly bad 
for small times as expected. The other parameters have (in 
both figures) been chosen as: 5r = 1, (5l = 2, eb. ~ sl ~ 0, 
£d = 1. 



Heisenberg picture) 



Sk = hkbk + hlbk'' 



dk = i{hkbk - Kbk^) (Al) 



a set of equations for the time-evolutions of (<t^), {(7^), 
{cr''}, {cr'^Sk}, {(T'^ak), (crysk), (cryak), {cr^Sk}, {cr^ak) 
in the Heisenberg picture, which is unfortunately not 
closed. However, we can achieve closure by assuming 
factorization of the expectation values and stationarity 
of the reservoir 



\ k' I 

( ^{sk'ak — ctfeSfc') ) 



2{a'^/y/^){sk') 

\hkf [1 + 2niuJk)] 
0, 



-2i\hi, 



(A2) 



which is consistent with the Born approximation. We can 
then analyze the steady state of the resulting equations 
and obtain the same results as discussed before, i.e., for 
pure dephasing {A = a^) we obtain (cr^oo) = (""^oo) = 
and {(T^oo) — (ctq) as discussed in subsection IIII C 21 and 
similarly for the dissipative case {A = a^) we obtain 

(o-'^oo) = (cr^oo) = and (cr^oo) = , which is 

consistent with Eqn. (|5T|) in subsection IIII C 31 



APPENDIX B: EXACT SOLUTION OF THE 
FANO- ANDERSON MODEL FOR LORENTZIAN 
TUNNELING RATES 



From the Hamiltonian we can calculate the time 
evolution of the fermionic operators in the Heisenberg 
picture (we use bold operator symbols to denote the 
Heisenberg picture and exploit the anticommutation re- 
lations throughout) 
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APPENDIX A: APPROXIMATE STEADY STATE 
OF THE SPIN-BOSON MODEL 

Starting from the Hamiltonian p4p we obtain with the 
abbreviations (we use bold symbols for operators in the 



k 

CkL = -iuJkhCkL + i^tkhd , 

CfcR = -lUJkRCkR + iMkRd , (Bl) 

which already forms a closed set of equations. These 
equations can be Laplace-transformed (where d{t) 
d{z) and d(t) — > —d + zd{z) and similarly for the other 
operators). 
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In Laplace space, we can eliminate Ckhiz) and CfcR(z) to solve the resulting equations for 

+ 00 



d{z) 



— 00 

(z + <5l + (z + <5r + )[d + ^\ + 



(z + ied)(z + <5l + ie-L){z + 5r + isR) + ^ [Tl5i^{z + (5r + ^er) + 5r(z + (5l + ieL)] ' 



(B2) 



where in the last line we have already assumed Lorentzian tunneling rates (|62p . The inverse Laplace transform 
(Bromwick integral [s^) can be performed by the theorem of residues d{t) = Ei R-6sJ(z)e"''^* , where Zi denote 

Z — Zi 

the poles of d{z). Denoting the roots of 

{z + iSd){z + (5l + isi^){z + (5r + is-R) + — [F^^^lI^ + (Jr + zsr) + ^\^i^{z + (5l + ^El)] = {z - zi){z - Z2){z - zaOBS) 

by zi, Z2, and Z3, respectively, we can easily calculate the residues (In case of degenerate roots, one may either use 
residue formulae for higher-order poles or simply analytic continuation of the solution for first order poles). For 
zi ^ Z2, zi ^ z^, and Z2 7^ Z3 we obtain the solution 



d{t) = 



(zi + (5l + zeL)(zi + (5r + zirR)e"'i* (z2 + 5l + ieL)(z2 + (^r + «eR)e^^* 



{zi - Z2)(zi - Z3) 

(z3 + (5l + jeL)(z3 + i5r + i£R)e'^3t 



(^2 - 2i)(z2 - Z3) 



n 



(2:3 - zi){zz - Z2) 

(zi + (5l + ieL)(2i + (5r + z£R)e^i* (z2 + (5l + «eL)(z2 + (^r + ieR)e''2* 



(21 - Z2)(zi - Z3)(zi + IWfcL) (2^2 - Zi)(z2 - Z3)(Z2 + iUkh) 



(z3 + (5L + ieL)(^3 + ^R + «£R)e'^='* , (-icj^L + + «£L)(-«t^fcL + <5r + i£R)e ^^^l* 



(Z3 - 2:i)(2:3 - ■Z2)(23 + i^k-L) {-ii^kh - Zi){-iuJkL - Z2){-iuJkL ~ Z3) 

(zi + (5l + ieL)(2i + ^R + i£R)e^i* (z2 + (5l + «eL)(^2 + <5r + ieR)e^=* 



+*A > , + 



(21 - Z2)(Z1 - Z3)(Z1 + iuJkR) (Z2 - Zl)(z2 - Z3)(z2 + ?WfcR) 



(z3 + (5l + ieL)(2;3 + + «eR)e^^* , (-icj^r + (5l + ieL)(-i'^feR + + «eR)e 



ifcLCfeL 



i:.RC,R. (B4) 



(^3 - Zl){z3 - Z2)(Z3 + ZWfcR) (-iWfcR - Zl)(-ZWfcR - Z2)(-iWfeR - Z3) 

With taking the initial conditions as (^cl,-^CkhJ = hk' hil^kh) and ^c[./RCfcR^ = 5uk' hi'^kY^) and {c\,^CkY^ = we 
obtain for nit) = {d^t)d{t)) the expression 

(zi + (5l + ieL)(2;i + + i£R)e^i* (z2 + (^l + i£L)(z2 + (^r + «eR)e^^* 



(Zl - Z2)(Z1 - Z3) 

(z3 + (Sl + i£L){z3 + Sk + «eR)e^^* 



(Z2 - Zi)(z2 - Z3) 



[zs - Zl)(z3 - Z2) 



no 



+ CXD 



+^ y [rL(c^)/LH + FRH/RH] X 

—00 

(zi + (5l + i£L)(zi + (5r + zeR)e^i* (z2 + (5l + jeL)(^2 + (^r + i£R)e^2* 



(21 - Z2)(zi - Z3)(zi + iuj) 



{Z2 - Zi){z2 - Z3)(Z2 + iuj) 



(z3 + (5l + ieL)(z3 + <5r + «£R)e^2* (-itj + (5l + «£l)(— it^ + 5r + i£R)e" 



(2^3 - zi)(z3 - Z2)(z3 + iuj) 



(—iuj — zi)(— — Z2){—iuJ — Z3) 



(B5) 



r 



In the large time- limit, this considerably simplifies (with Fl(i-li) and Fr(ci;) and by setting A — > 1 we explicitly 
using that dizi < 0). Conventionally, A^ is absorbed in 
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recover the well-known steady-state results in the liter- 
ature (compare e.g., Eqns (12.27) with using Lorentzian 
tunneling rates of form dM]), Eqns. (12.30), and (12.31) 
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